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THE  VIEW,  OPINIONS,  AND/OR  FINDINGS  CONTAINED  IN  THIS  REPORT  ARE 
THOSE  OF  THE  AUTHORS)  AND  SHOULD  NOT  BE  CONSTRUED  AS  AN  OFFICIAL 
DEPARTMENT  OF  THE  ARMY  POSITION,  POLICY,  OR  DECISION,  UNLESS  SO 
DESIGNATED  BY  OTHER  DOCUMENTATION. 


Abstract 

A  model  of  eolian  sediment  transport  has  been  constructed,  a  special  case 
of  which  is  that  corresponding  to  sand-sized  mineral  grains  subjected  to 
moderate  winds:  saltation.  The  model  consists  of  four  compartments 
corresponding  to  (1)  aerodynamic  entrainment,  (2)  grain  trajectories,  (3) 
grain-bed  impacts,  and  (4)  momentum  extraction  from  the  wind.  Each 
sub-model  encapsulates  the  physics  of  the  process,  and  is  constrained 
where  necessary  by  experimental  data.  When  combined,  the  full  model 
allows  simulation  of  eolian  saltation  from  inception  by  aerodynamic 
entrainment  to  steady  state. 

Many  observed  characteristics  of  natural  saltation  systems  are  reproduced 
by  the  simulations.  Steady  state  mass  flux  and  concentration  profiles  all 
display  rapid  decay  with  height  above  the  bed,  representing  the 
preponderance  of  short,  low-energy  trajectories  in  the  saltation  population. 
Yet  the  role  of  less  abundant,  longer,  higher-energy  trajectories  is  a  strong 
one:  at  steady  state  the  entire  population  of  saltating  grains  is  controlled  by 
bed  impacts  rather  than  aerodynamic  entrainment.  Because  the  nature  of 
the  grain  splash  process  is  such  that  high-energy  impacts  are  much  more 
efficient  at  ejecting  other  grains  from  the  bed,  the  response  time  of  the 
system  to  changes  in  wind  velocity  is  determined  by  the  hop  time  of  these 
long  trajectories.  Several  hop  times,  or  roughly  1-2  seconds,  are  required. 

Varying  wind  velocity  among  the  simulation  runs  allows  mapping  of  the 
relation  between  steady  state  mass  flux  and  wind  velocity  —  the  mass  flux 
"law"  --  which  may  be  expressed  as  a  power  law  of  the  excess  shear 
velocity.  The  hysteresis  that  led  Bagnold  to  define  both  fluid  and  impact 
thresholds  for  saltation  is  apparent,  reinforcing  our  conclusion  that  it  is  the 
impacts  of  saltating  grains  that  supports  the  large  population  of  saltating 
grains  at  steady  state. 

These  characteristics  point  to  a  serious  need  to  collect  high  frequency  near¬ 
bed  wind  velocity  data  simultaneously  with  mass  flux  measurements  in 
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natural  environments.  Ideally  this  data  should  be  collected  at  a  rate  of 
several  Hz,  in  order  to  capture  the  relevant  spikes  in  wind  velocity  that 
contribute  most  strongly  to  the  total  flux  over  any  period  of  time  long 
relative  to  the  response  time  of  the  saltation  system.  Sample  calculated 
time  series  of  both  wind  speed  (constrained  by  empirically  gathered  wind 
spectra)  and  the  resulting  saltation  flux  (taking  into  account  both  the 
response  time  and  the  mass  flux  law)  demonstrate  the  importance  of  high 
frequency  variations  in  wind  speed. 

Although  such  models  lack  the  detail  to  treat  the  true  complexity  of  the 
natural  world  they  serve  both  to  develop  our  intuition  about  the  natural 
setting,  and  to  guide  future  experimental  efforts. 
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Introduction 


Considerable  progress  has  been  made  in  the  understanding  of  the  process 
by  which  sand  sized  grains  are  transported  by  the  earth's  atmosphere. 
Although  ideas  initially  espoused  by  early  workers,  especially  Bagnold 
(1941)  and  Chepil  (1945),  have  enjoyed  broad  support  over  the  decades 
through  empirical  findings,  only  recently  has  the  theory  of  sand  transport 
developed  to  provide  the  firm  underpinnings  to  such  empirical  work. 
Theoretical  work  has  now  reached  an  approximate  par  with  the  empirical, 
and  new  questions  raised  by  this  theoretical  work  will  require  new  field 
and  wind  tunnel  efforts  in  the  future. 

The  broad  motivations  for  the  study  of  the  simple  system  of  sand  plus  air 
range  widely,  but  may  be  classed  in  four  camps:  environmental,  geological, 
planetary  and  physical.  The  first  motivated  the  work  of  Chepil  and  others, 
who  lived  in  and  through  the  Dust  Bowl  years  of  the  US.  Research  founded 
by  Chepil  has  been  carried  out  since  in  extensive  field  and  wind  tunnel 

work  centered  at  the  Kansas  research  station,  with  the  major  emphasis 
being  placed  on  the  effects  of  wind  on  crops  of  the  Great  Plains  and  how  to 
minimize  these  effects. 

The  geological  motivations  for  performance  of  research  on  eolian  sand 
systems  include  the  need  to  understand  the  mechanisms  by  which  eolian 

sandstone  bodies  of  high  porosity  are  produced,  and  how  the  grain  packing 
geometry  and  other  intrinsic  variables  of  importance  to  the  oil  industry 

might  differ  according  to  the  important  variables  in  the  problem:  wind 
velocity,  grain  size  distribution,  etc. 

The  strong  motivation  to  understand  the  settings  into  which  we  were 

sending  planetary  landing  craft,  made  the  1970's  and  early  1 980’s  ripe  for 
a  resurgence  of  work  on  both  the  theoretical  and  empirical  sides  of  the 
problem.  Wind  tunnels  were  constructed  capable  of  very  high  and  very 
low  atmospheric  pressures,  in  which  both  threshold  velocities  and  mass 
flux  relations  were  measured;  high-speed  filming  techniques  were  used  to 
track  individual  trajectories  in  saltation;  interest  was  renewed  in  terrestrial 
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settings  in  which  analogues  to  Martian  terrains  were  suggested  by  surface 
features;  and  theoretical  work  was  renewed  on  both  grain  trajectories  and 
threshold  velocities  for  a  wide  range  of  atmospheric  conditions. 

The  fourth  motivation  may  be  termed  physical  in  that  the  research  has 
been  on  the  physics  of  the  problem,  rather  than  on  its  planetary,  geological, 
or  environmental  applications.  The  several  components  of  the  problem 
have  been  taken  apart  and  studied  in  detail:  e.g.,  the  grain  trajectory,  the 
grain-bed  impact,  the  wind  velocity  feedback,  and  the  aerodynamic 
threshold.  Only  recently  have  computational  capabilities  allowed  putting 
all  of  these  components  back  together  again.  In  addition,  interest  was 
spurred  by  the  degree  to  which  the  sand-air  system  was  self- 
organizational,  resulting  in  sand  ripples  of  a  well-defined  wavelength. 

Below  we  give  an  overview  of  the  current  state  of  knowledge  of  the 
saltation  problem,  emphasizing  results  obtained  as  a  result  of  Army 
Research  Office  support. 


The  physical  picture 

It  is  useful  to  picture  eolian  saltation  in  its  simplest  manifestation  as  being 
comprised  of  four  linked  processes:  aerodynamic  entrainment,  grain 
trajectories,  grain-bed  impacts  (the  granular  splash  problem),  and  the  wind 
field  modification.  Clearly,  no  sand  is  transported  at  the  very  lowest  wind 
velocities.  Sand  transport  begins  with  increasing  wind  velocity  as  the  most 
easily  moved  grains  are  displaced.  The  aerodynamics  of  the  threshold 
condition  must  therefore  be  known,  and  quantified  appropriately.  Once 
airborne,  grains  travel  trajectories  that  are  determined  by  the  balance 
between  grain  weight,  drag  and  lift  forces  on  the  grain.  The  smaller  the 
grains,  the  more  likely  they  are  to  respond  effectively  to  the  turbulence  of 
the  boundary  layer.  We  will  consider  in  this  paper  only  those  grains  whose 
size  is  large  enough  to  filter  out  short  time  scale  variations  in  wind  velocity, 
meaning  that  they  respond  effectively  only  to  the  mean  wind  velocity 
profile:  we  therefore  limit  the  discussion  to  saltation,  rather  than 

suspension  (see  Anderson,  1987).  Importantly,  this  reduces  the  grain 
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trajectory  problem  to  a  deterministic  one:  once  the  liftoff  conditions  are 
specified,  the  entire  grain  trajectory  is  specified,  including  among  other 
things  the  hop  height,  length,  and  hop  time.  Due  to  gravitational 
acceleration,  all  grains  that  are  ejected  from  the  bed  ultimately  return  to 
the  bed,  upon  which  they  interact  energetically  with  a  collection  of  other 
nearby  particles.  We  must  know  in  a  quantitative  sense  the  nature  of  this 
splash  process,  including  what  happens  to  the  impacting  grain,  and  to 
nearby  grains.  Because  of  the  complicated  microtopography  of  even  a  "flat" 
bed  and  because  of  the  complexity  of  the  packing  geometry,  this  splash 
process  is  stochastic:  we  may  know  the  results  of  any  particular  impact, 
characterized  by  the  size  of  the  impacting  grain  and  its  angle  and  speed, 
and  the  nature  of  the  local  bed  (grain  size  distribution,  angle  with  respect 
to  horizontal),  in  only  a  statistical  sense.  All  of  the  above-mentioned 
aspects  of  the  saltation  problem  are  dependent  only  upon  the  behavior  of 
individual  grains.  The  fourth  facet  to  the  saltation  problem,  that  of  the 
modification  of  the  wind  velocity,  reflects  the  extraction  of  momentum 
from  the  wind  by  ail  of  the  particles  in  the  air,  and  is  an  integral  quantity 
that  requires  knowing  the  behavior  of  all  of  the  airborne  particles.  In  turn 
this  modification  will  affect  all  of  the  other  three  processes,  and  may 
therefore  be  thought  of  as  a  negative  feedback  mechanism:  aerodynamic 
entrainment  will  be  curtailed  to  some  extent,  grain  trajectories  will  be 
shortened  by  the  fact  that  the  entire  wind  velocity  profile  is  slowed,  and 
because  the  grain  trajectories  are  slowed  the  impact  velocities  are 
diminished.  Calculations  such  as  we  present  here,  which  extend  the  work 
presented  elsewhere  (Anderson  and  Haff,  1988)  allow  determination  of  the 
relative  importance  of  these  two  feedback  effects. 

We  will  first  introduce  each  of  the  four  processes  mentioned,  with  attention 
being  paid  in  inverse  proportion  to  attention  in  the  present  literature.  We 
then  present  the  results  of  calculations  performed  with  all  processes  linked. 

Aerodynamic  entrainment 

Grains  may  be  entrained  by  direct  aerodynamic  forces  of  lift  and  drag. 
Aerodynamic  entrainment  is  required  for  the  initiation  of  saltation  in  the 
absence  of  external  mechanical  disturbances.  As  the  simplest  possible 
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model  for  aerodynamic  entrainment,  we  take  the  number  of  entrained 
grains  per  unit  area  of  bed  per  unit  time  (the  aerodynamic  ejection  rate)  to 
be  proportional  to  the  excess  shear  stress 


Na  =  a  (xa-xc) 


where  xa  is  the  short  term  mean  shear  stress  at  the  bed,  xc  is  the  critical 

fluid  shear  stress  for  entrainment,  for  which  an  expression  has  been 
determined  from  numerous  wind  tunnel  experiments,  and  a  is  a  constant. 
All  such  aerodynamically  entrained  grains  are  assumed  to  leave  the  bed 
with  a  velocity  corresponding  to  that  needed  to  reach  a  height  of  about  one 
grain  diameter,  D,  above  the  bed.  This  is  a  conservative  assumption  that 
effectively  limits  the  role  of  the  aerodynamic  grains  in  the  eolian  saltation 
process. 

As  will  be  demonstrated  below,  this  somewhat  ad  hoc  entrainment  rule 
does  not  play  a  critical  role  in  a  well-developed  saltation  curtain. 


Grain  trajectories 

Once  launched  into  the  airstream  by  either  ballistic  impact  forces  or 
aerodynamic  forces,  grains  are  acted  upon  by  the  surface  tractions  of  drag 
and  lift  forces  imposed  by  the  wind,  and  by  the  body  force  of  the  grain's 
weight.  Many  studies  have  addressed  the  details  of  the  forces  acting  upon 
airborne  grains,  beginning  with  Bagnold  (1941).  They  may  be  simply 
stated  as  follows.  Consider  a  coordinate  system  in  which  positive  x  is  in  the 
downwind  direction  and  positive  z  is  upward.  The  particle  weight  is  Fw, 

written 


Fw  -  Mp  g 

where  Mp  is  the  particle  mass  and  g  the  acceleration  due  to  gravity.  The 

drag  force,  Fd,  operates  in  the  direction  of  the  relative  velocity  of  the  grain 
and  the  local  wind  velocity,  Urel,  and  may  be  written 
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Fd  =  ^PaCdAIUre|l  Ure| 

where  p a  is  the  density  of  air  (-1.22kg/m3),  A  is  the  cross-sectional  area  of 
the  grain  presented  to  the  flow  (=7tD2/4  for  spheres,  where  D  is  grain 
diameter),  and  C<j  is  the  drag  coefficient,  which  is  dependent  upon  the  local 
Reynold's  number,  Re=UrelD/v,  where  v  is  the  kinematic  viscosity  of  air 
(=1.5xlO-5  m2/s  at  20°C). 

As  the  shear-imposed  lift  force  on  a  spherical  grain  is  essentially  negligible 
above  a  few  grain  diameters  above  the  bed  (Anderson  and  Hallet,  1986), 
we  may  neglect  it  in  the  first  order  sketch  of  the  trajectory  problem 
presented  here.  More  elaborate  and  rigorous  treatments  will  incorporate 
this  effect  as  well  as  others,  among  them  the  Magnus  effect  (White  and 
Schulz,  1977)  and  the  virtual  mass  effect.  The  essence  of  the  trajectory 
problem  is,  however,  contained  in  the  weight  and  drag  forces.  The 
downstream  and  vertical  components  of  these  two  forces  may  be  added  to 
yield  the  resultant  forces  in  the  x  and  z  directions,  from  which 
accelerations,  velocities  and  displacements  may  be  calculated  for  a  given 
time  step,  dt,  from  Newton's  second  law: 

a=djl  =  Fw+Fd  u=d*  =  dlldt;  X  =  dx-dt 

dt  Mp  dt  dt  dt 

Given  the  initial  conditions  of  liftoff  velocity  and  speed,  these  equations  can 
be  integrated  forward  in  time  until  the  particle  reimpacts  the  bed.  We  find 
that  very  accurate  calculations  of  trajectories  may  be  accomplished  with 
100  time  steps  per  trajectory,  requiring  that  the  time  step,  be  chosen  such 
that  dt=xhop/100,  where  Thop  is  the  expected  hop  time  in  the  absence  of 
vertical  drag,  or  2w0/g,  w0  being  the  initial  vertical  velocity  of  the  grain. 

[It  is  the  happy  case  in  the  eolian  trajectory  problem  that  the  accelerations 
in  the  vertical  are  due  primarily  not  to  fluid  drag  but  to  gravity,  meaning 
that  hop  times  and  hop  heights  are  well  approximated  (less  about  10-20%) 
by  their  values  in  a  vacuum  (see  Anderson  and  Hallet,  1986).] 
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The  splash  process 


Our  knowledge  of  the  splash  process  has  been  augmented  greatly  in  the  last 
decade  by  both  physical  and  numerical  experiments.  The  sophistication  of 
these  experiments  has  advanced  considerably,  and  has  resulted  in  placing 
us  in  a  position  to  characterize  in  considerable  detail  the  complexity  of  the 
process.  Saltation  is  a  stochastic  process  that  must  take  into  account  the 
reaction  not  only  of  the  grain  impacting  the  bed,  but  the  resulting 
rearrangement  of  the  grains  in  the  bed,  including  the  possibility  of  ejecting 
some  of  these  grains.  The  impact  process  is  expected  to  depend  upon  the 
mass  and  velocity  of  the  impacting  grain,  including  its  speed  and  approach 
angle  with  respect  to  the  local  bed  slope,  and  the  masses  and  elastic 
parameters  describing  the  grains  comprising  the  bed  near  the  impact  point. 

Early  numerical  experiments  were  crude  in  the  sense  that  they  took  into 
account  only  the  geometrical  complexity  of  the  problem,  while  ignoring  the 
dynamical  nature  of  the  problem.  Grains  in  the  bed  were  not  allowed  to  be 
ejected,  and  the  splash  process  was  reduced  to  finding  the  distribution  of 
rebound  angles  to  be  expected  from  beds  of  particular  micro-topographies. 
Simple  geometrical  models  such  as  these  predict  that  if  the  bed  were 
perfectly  rigid  the  rebounding  particle  would  most  likely  emerge  from  the 
collision  at  roughly  50°  from  the  horizontal,  having  impacted  the  bed  at 
angles  typical  of  eolian  saltation,  i.e.  10-15°.  This  matches  closely  the 
measured  averages  of  saltation  angles  in  full  saltation  experiments,  as  first 
reported  by  White  and  Schulz  (1977),  and  more  recently  in  numerous 
publications  by  Willetts  and  Rice  (e.g.  1986a,b  ,1988,  1989).  This  was  also 
the  first  work  to  question  seriously  the  conclusions  of  Bagnold  (1941)  that 
the  preponderance  of  grains  were  ejected  from  the  bed  vertically, 
conclusions  that  had  played  an  important  role  in  shaping  later  theoretical 
work  (e.g.  Owen,  1964). 

The  utilization  of  high-speed  films  by  Willetts  and  Rice  at  the  University  of 
Aberdeen  has  allowed  the  characterization  of  the  splash  problem  during 
saltation  in  wind  tunnel  settings.  By  elevating  a  small  strip  of  the  bed  out 
of  the  high-concentration  portion  of  the  saltation  layer,  they  were  able  to 
film  individual  impacts  from  which  a  statistical  description  of  the  impact 
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process  could  be  gathered.  Typical  saltation  impacts  strike  the  bed  at  angles 
of  10-15°  from  the  horizontal.  The  impact  results  in  two  phenomena:  the 
rebound  of  the  impacting  grain,  and  the  emergence  of  several  grains  from 
among  the  bed  grains.  The  rebound  speed  is  reduced  by  roughly  half  from 
the  pre-impact  value,  and  the  mean  rebind  angle  is  roughly  50°  from  the 
vertical.  The  number  of  splashed  grains  rises  roughly  linearly  with  the 
impact  speed,  being  of  the  order  of  3  grains  per  impact  for  impact  speeds 
on  the  high  end  of  those  expected  in  eolian  saltation.  The  ejection  speeds  of 
the  splashed  grains  are  of  the  order  of  10  percent  of  the  impact  speed,  and 
the  ejection  angles  are  canted  more  steeply  toward  the  vertical  than  the 
rebound,  within  a  broad  range. 

These  measurements  have  been  used  to  check  the  results  of  numerical 
simulations  that  attempt  to  take  into  account  all  the  forces  between  all  the 
grains  in  the  problem.  This  was  first  addressed  in  detail  by  Werner  (e.g., 
1987;  1990)  in  his  dynamical  simulations  of  the  impact  process.  Working 
with  well-packed  two  dimensional  beds,  impactor  grains  could  be  fired  into 
the  beds  at  any  chosen  angle  and  speed.  The  statistics  from  many  such 
impacts  were  collected  to  describe  the  outcome  from  any  set  of  impact 
parameters.  Werner  reported  distributions  of  rebound  and  ejection  angles 
that  paint  a  picture  similar  to  that  obtained  from  the  high  speed  films. 
Most  of  the  grains  ejected  from  the  bed  were  found  to  emerge  from  within 
the  top  one  to  two  grain  layers,  and  displayed  a  high  ejection  angle. 
Werner  noted  that  this  was  probably  due  to  the  fact  that  this  was  the 
direction  of  the  free  face,  and  therefore  the  direction  in  which  there  was 
least  resistance  to  motion.  He  further  noted  that  when  the  bed  was  given 
some  micro-relief,  that  grains  were  preferentially  ejected  from  micro¬ 
discontinuities  he  termed  steps.  This  allowed  a  greater  range  of  possible 
ejection  angles,  and  therefore  resulted  in  a  greater  spread  of  the  resulting 
ejection  angle  distribution.  Simultaneously,  the  Caltech  group  (e.g.,  Mitha  et 
al.,  1986)  developed  a  means  of  producing  and  documenting  the  results  of 
individual  grain  impacts  using  real  particles.  These  experiments  started 
using  spherical  steel  pellets  (bb's)  which  were  fired  at  a  given  angle  into  a 
wide  bed  of  similarly  sized  pellets.  Subsequent  development  of  a  "sand 
gun"  allowed  similar  experiments  using  individual  sand  grains  fired  into 
sand  beds  (Werner,  1987). 
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Subsequent  to  Werner's  simulations,  Haff  and  Anderson  (1988)  reported 

further  dynamical  simulations  that  attempt  to  map  out  the  full  splash 

function.  The  well-packed  nature  of  the  simulated  bed  is  abandoned  in 
these  calculations.  Several  two-dimensional  simulated  beds  are  created  by 
dropping  grains  from  a  given  height  with  a  random  distribution  of  initial 
velocities.  The  particles  are  allowed  to  settle  in  a  "box"  that  is  constructed 
with  wrap-around  boundary  conditions  that  effectively  extend  the  bed 
laterally.  The  resulting  beds  manifest  more  realistic  micro-topography  that 
importantly  affect  the  splash  outcome.  These  beds  are  then  each  subjected 
to  impacts  with  identical  impact  parameters  of  speed  and  angle,  many 
impacts  being  performed  on  each  of  three  beds,  each  impact  impinging  on 
the  bed  at  a  different  point  in  order  to  collect  sufficient  statistics.  The 
resulting  splash  descriptions  therefore  have  embedded  in  them  the 

stochastic  nature  of  the  bed,  as  well  as  the  possibility  of  hitting  any 

particular  grain  in  the  bed  at  any  of  many  different  points. 

Impact  model 

The  splash  function  used  in  these  calculations  is  generated  from  direct 
computer  simulation  of  the  motion  of  the  impacting  and  impacted  grains. 
The  simulation  is  performed  in  two-dimensions,  with  each  grain  assumed  to 
be  circular.  The  position  and  velocity  of  each  particle  in  the  simulation  is 
followed  explicitly  through  time  via  integration  of  the  Newtonian  equations 
of  motion 

Fj  =  Mj  a; 

where  Fj  is  the  force  applied  to  particle  i,  Mj  is  the  particle  mass,  and  aj  is 
its  acceleration.  Since  friction  between  particles  is  included  in  the  model, 
torques  may  be  applied  leading  to  particle  rotation.  To  the  three  degrees  of 
freedom  attaching  to  the  motion  of  a  two-dimensional  extended  particle, 
namely  translational  motion  in  the  x  and  z  directions,  plus  angular  motion 
about  the  center  of  mass,  there  correspond  the  two  components  of  Newton's 
second  law,  Fx=Max  and  Fz=Maz  plus  a  torque  equation, 
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where  x  is  the  applied  torque,  I  the  moment  of  inertia  about  the  center  of 
mass,  and  a  the  particle's  angular  acceleration.  Following  the  motion  of  each 
particle  therefore  requires  the  integration  of  3  second  order  or  6  first  order 
differential  equations,  for  a  total  of  6N  equations  for  a  model  impact  system 
with  N  grains.  Typical  impact  scenarios  have  involved  80  to  100  grains, 
corresponding  to  ~  500  -  600  simultaneous  equations. 

There  are  two  basic  types  of  forces  which  we  must  consider,  body  forces 
and  contact  forces.  Contact  forces  can  be  categorized  as  normal  (or 
compressive),  and  tangential  (or  frictional).  The  relevant  body  force  here  is 
just  the  weight  of  the  grain,  Mg  (we  drop  the  subscript  i  for  now). 

The  normal  force  is  modeled  as  a  very  stiff  linear  spring.  This  force  is  not 
activated  until  two  particles  "contact"  one  another,  as  determined  by  the 
spacing  of  their  centers  becoming  less  than  the  sum  of  their  radii.  The 
stiffness  of  the  spring  is  chosen  great  enough  that  particle  "overlap"  is 
always  less  than  a  predetermined  fraction  of  the  particle  diameter,  a 
consideration  which  depends  upon  the  maximum  relative  particle  velocities 
expected  to  occur  in  the  simulation.  In  the  saltation  problem,  this  is  always 
determined  by  the  speed  of  the  impacting  grain.  Higher  speeds  require 
stiffer  springs.  Since  the  natural  period  of  a  linear  spring  decreases  as  the 
spring  constant  increases, 

T  =  VM7¥ 

an  increase  in  k  implies  an  increase  in  the  number  of  incremental  time 
steps  At  which  are  needed  to  effect  the  actual  integration  of  the  equations 
of  motion,  (since  we  must  have  A  t « T  ).  This  is  one  limitation  on  the 
number  of  independent  splash  simulations  which  can  in  practice  be  run,  i.e., 
T  must  be  small  compared  to  the  time  required  for  the  bed  to  respond 
sufficiently  that  final  grain  trajectories  can  be  determined.  And  A  t  is 
generally  taken  to  be  at  least  an  order  of  magnitude  smaller  in  order  to 
ensure  adequate  accuracy  of  integration.  (The  integration  itself  is  carried 
out  using  a  predictor-corrector  scheme  for  which  preset  integration 
tolerances  can  be  specified.) 
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This  calculation  is  made  slightly  more  complicated  by  the  fact  that  the  grain 
collision  is  be  damped,  i.e.,  the  grains  are  slightly  inelastic.  Damping  is 
modeled  by  a  velocity- dependent  force  which  is  proportional  to  the 
instantaneous  normal  relative  particle  velocity  and  which  always  opposes 
relative  particle  motion.  In  such  a  damped  oscillator,  there  are  two  time 
scales,  one  is  the  effective  or  modified  spring  frequency,  which  is  less  than 
the  natural  (undamped)  frequency,  and  the  other  is  the  damping  time  scale 
which  gives  the  e-folding  time  for  degradation  of  the  oscillator  amplitude. 
For  an  accurate  integration  of  the  equations  of  motion,  the  integration  time 
step  A  t « T  must  be  less  than  both  of  the  system  physical  time  scales. 

For  such  an  oscillator  there  are  three  characteristic  response  modes, 
namely,  underdamping,  critical  damping  and  overdamping,  depending  upon 
whether  the  damping  time-scale  is  greater  than,  equal  to,  or  less  than  the 
natural  frequency  of  the  system.  Critical  and  overdamped  behavior 
correspond  to  collisions  in  which  the  two  colliding  partners  lose  all  their 
normal  velocity  in  a  head-on  collision  (i.e.,  the  coefficient  of  restitution,  e, 
given  by  the  ratio  of  outgoing  to  incoming  (normal)  relative  velocity,  is 
zero).  For  the  collision  of  resilient  quartz  grains  the  coefficient  of  restitution 
is  greater  than  zero  (it  would  be  unity  for  a  perfectly  elastic  collision)  and 
consequently  the  normal  forces  in  the  impact  simulation  are  modeled  by 
underdamped  springs.  The  choice  of  spring  constant  and  damping  strength 
is  therefore  made  to  satisfy  the  overlap  requirements,  and  to  produce  a 
desired  coefficient  of  restitution.  The  precise  value  of  e  to  be  used  in  a 
given  simulation  is  somewhat  problematical,  because,  at  the  collision  speeds 
characterizing  typical  impact  events  (a  few  ms’1),  the  actual  grain  contact  is 
not  at  all  elastic,  but  rather  is  characterized  by  plastic  deformation,  spalling 
and  chipping  off  of  small  pieces  of  the  parent  grains.  Simulation  runs  on 
identical  systems  with  reasonable  variation  of  e  show  that  the  main  details 
of  the  splash  function  are  not  very  sensitive  to  the  exact  value  of  e.  In  our 
previous  simulations  we  have  taken  e  =  0.7. 

Besides  the  normal  compressive  and  the  corresponding  damping  forces, 
which  act  in  a  direction  along  the  line  of  centers  of  two  contacting  particles, 
there  is  also  a  tangential  force  due  to  friction  at  the  point  of  contact.  Contact 
friction  is  modeled  in  one  of  two  different  ways,  depending  upon  whether. 
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at  any  instant,  there  is  slippage  at  the  contact  or  not.  If  the  contact  is 
slipping,  i.e.,  if  there  is  relative  tangential  motion  of  the  material  just  on 
either  side  of  the  contact  point,  we  invoke  a  rule  corresponding  to  Coulomb 
friction,  i.e.,  that  the  tangential  frictional  force  is  related  to  the 
(instantaneous)  normal  force  Fn  via  a  constant  of  proportionality,  the 
friction  coefficient  4,  If  the  contact  is  not  slipping,  then  we  imagine  that  a 
tangential,  damped,  linear  spring  acts  to  oppose  any  slippage.  Particles  may 
still  roll  freely  because  the  friction  spring  is  a  function  of  the  total  arc 
length  separating  the  anchoring  points  of  the  spring  minus  the  amount  of 
rolling.  If  slippage  starts  to  occur,  the  very  stiff  tangential  spring 
effectively  limits  rim  displacements  unless  or  until  the  spring  force  exceeds 
4  Fn  .  If  this  occurs,  then  the  tangential  spring  is  severed,  and  the  maximum 
tangential  force  remains  at  4  Fn.  The  tangential  spring  is  damped  in  order  to 
suppress  low-amplitude,  high-frequency  (unphysical)  tangential  oscillations 
of  the  contacting  particles. 

The  choice  of  a  friction  coefficient  is  somewhat  problematical  also,  because, 
while  perfectly  smooth  quartz  grains  exhibit  low  friction,  real  sand  grains 
are  pitted,  faceted  and  generally  irregular  in  shape.  Thus,  in  nature,  torques 
can  be  exerted  across  a  particle  contact  even  in  the  absence  of  true  friction 
(Clayhold  effect).  The  friction  coefficient  4,  for  the  disk  model  described 
here,  is  clearly  an  effective  friction  coefficient  which  incorporates  not  only 
the  effects  of  true  surface  friction,  but  also  effects  of  particle  shape  and 
surface  morphology.  In  our  saltation  calculations,  we  have  used  a  generic 
value  of  0.5  for  the  friction  coefficient.  As  with  the  coefficient  of  restitution, 
the  splash  function  is  only  modestly  sensitive  to  4,  as  long  as  it  is  not  too 
small. 

Each  saltation  impact  calculation  is  carried  forward  in  time  far  enough  that 
the  nature  of  the  bed-response  is  well-established.  Depending  upon  bed 
conditions  this  might  be  10  -20  T.  If  the  simulation  is  followed  too  far  in 
time,  complications  may  arise  from  the  presence  of  stress  waves  reflected 
off  the  substratum  upon  which  the  simulated  collection  of  bed  grains  rests. 
These  waves,  which  are  artifacts  of  the  calculation,  can  actually  throw 
grains  off  the  surface  and  lead  to  an  overestimation  of  the  number  of 
reptating  grains.  This  is  one  of  several  computational  issues  which  arise  as 


a  result  of  the  necessarily  small  number  of  grains  which  can  be  handled  in 
each  impact  event. 

In  our  simulation  studies  the  "target”  bed  was  prepared  by  dropping 
between  80  and  100  grains  onto  a  fixed  and  immobile  surface.  The  grains 
were  allowed  to  fall  under  the  influence  of  gravity  and  to  rebound  back 
again,  colliding  with  one  another  however  they  might,  until  friction  and  the 
inelasticity  of  collisions  reduced  the  kinetic  energy  of  the  grain  mass  to 
essentially  zero.  Because  of  the  relatively  small  number  of  grains,  periodic 
boundary  conditions  were  applied  at  the  edges  of  the  region  of  calculation. 
At  these  lateral  boundaries,  particles  and  their  contacts  are  "wrapped" 
around  to  the  other  side  of  the  cell,  so  that  if  a  particle  moves  too  far  to  the 
right  and  crosses  the  right-hand  dashed  line,  it  is  inserted  again  on  the  left, 
at  precisely  the  same  elevation  it  had  on  the  right.  While  periodic  boundary 
conditions  can  introduce  spurious  behavior  in  certain  situations  (especially 
in  particle  flow  studies),  they  seem  to  be  relatively  innocuous  here  because 
the  calculation  of  every  impact  is  wisely  terminated  before  the  disturbance 
initiated  by  the  impact  can  propagate  out  one  side  of  the  test  region  and 
back  in  the  other. 

Once  a  bed  has  been  constructed  in  this  way,  we  impinge  an  identical 
projectile  of  given  velocity  and  impact  angle  upon  it.  Ten  impact  positions 
were  chosen  for  each  velocity  and  impact  angle  combination,  corresponding 
to  equally  spaced  impact  points  for  a  reference  flat  bed.  Since  the  bed  has 
its  own  unique  topography,  the  actual  impact  points  were  not  uniformly 

distributed.  Thus  the  "upwind"  side  of  any  clumping  of  grains  at  the  surface 
will  tend  to  receive  somewhat  more  impacts  than  the  "downwind"  side  of 
such  a  clump,  as  would  be  the  case  in  nature.  (This  is  the  physical  origin  of 
wind-ripple  instabilties.) 

One  technical  issue  of  importance  in  grain  simulation  revolves  around 
contact  detection.  Two  grains  are  in  contact  whenever  the  separation 

between  their  centers  of  mass  is  less  than  the  sum  of  their  radii.  Since  all 

particles  in  the  simulation  may  potentially  move,  we  must  have  a  method 

for  periodically  detecting  the  generation  of  new  contacts  and  the  breaking 
of  old  ones.  One  way  to  do  this  is,  at  every  integration  step,  to  check  the 
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separation  of  the  centers  of  mass  of  every  pair  of  particles.  If  there  are  N 
particles  in  the  system,  then  there  are  N(N-l)/2  pairs,  so  the  number  of 
pairs  to  be  checked  increases  like  N2  for  large  N.  Implementing  contact 
detection  this  way  leads  to  a  rapid  saturation  of  computing  power  when  the 
number  of  particles  becomes  large.  Various  techniques  exist  to  cut  down 
the  fraction  of  computer  time  spent  in  contact  detection.  Most  of  them 
involve  some  sort  of  fine-  or  cross-graining  and  take  advantage  of  the  fact 
that  a  grain's  neighbors  (contacting  or  not)  at  time  t  will  still  be  the  grain’s 
(only)  neighbors  at  time  t  +  A  tcheck  ,  where  tcheck  is  large  compared  with  the 
integration  time  step  A  t,  but  small  compared  with  the  total  elapsed  time  in 
the  simulated  system,  ttotai  .Thus,  during  an  interval  of  time  A  tcheck  it  is 
necessary  to  check  only  n « N  contacts  within  a  small  neighborhood  of  each 
particle.  Only  occasionally  (at  tcheck  intervals)  need  all  pairs  be  checked. 

The  eolian  saltation  impact  program  implements  a  variation  of  this  method, 
in  that  each  particle  carries  attached  to  it  a  "neighbor  list",  which  contains 
an  identifying  particle  number  for  every  (other)  particle  that  was  within  a 
specified  radius  at  time  t.  At  each  At,  the  collision  detection  code  looks  only 
at  the  neighbor  list.  The  list  is  updated  from  time  to  time,  using  a  rule 
based  upon  the  speed  of  the  fastest  particle  in  the  system,  as  old  neighbors 
move  away  and  new  ones  move  in. 

Numerical  simulations  of  single  grain  impacts  into  a  granular  bed  were 
performed  with  the  aim  of  providing  a  quantitative  picture  of  the  splash 
process.  The  impact  simulations  were  performed  in  two  dimensions  with 
identical  slightly  inelastic  1  mm  diameter  grains  (disks)  whose  interactions 
were  characterized  by  a  coefficient  of  restitution  equal  to  0.7,  an  inter¬ 
grain  friction  coefficient  equal  to  0.5,  and  an  elastic  modulus  equal  to  10^ 
dy  cm'  1 .  Beds  were  generated  by  dropping  87  particles  with  random 
initial  velocities  into  a  box  with  periodic  boundaries,  its  width  chosen  such 
that  the  resulting  bed  was  about  10  grains  deep.  Three  such  granular  beds, 
distinguished  by  different  packing  configurations,  were  each  impacted  20 
times  to  develop  the  splash  statistics  for  each  impact  angle  and  speed. 
These  computations  soften  assumptions  made  in  earlier  numerical  work  on 
impact  dynamics  (Werner  1987). 
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The  picture  of  the  grain  splash  resulting  from  our  numerical 
simulations  is  qualitatively  similar  to  that  derived  from  earlier  numerical 
(Werner  1987,  1990)  and  physical  experiments  (Mitha  1986;  Willetts  and 
Rice  1989)  and  fills  in  gaps  in  these  latter  results,  especially  in  the  low- 
impact  velocity  range.  In  our  simulations  the  impacting  grains  were  found 
to  rebound  from  the  surface  with  nearly  unit  (0.95)  probability,  a  result 
postulated  earlier  by  Rumpel  (1985).  For  the  chosen  grain  parameters,  the 
mean  rebound  speed  is  approximately  50-60%  of  the  impact  speed,  and  the 
mean  rebound  angle  is  30-40°  from  the  horizontal.  The  impact  forces 
result  in  ejection  of  a  number  of  grains  from  within  a  few  grain  diameters 
of  the  impact  site.  The  mean  number  of  grains  ejected  increases  roughly 
linearly  with  impact  speed,  and  compares  well  with  data  obtained  from 
physical  impact  experiments  with  coarse  sand.  The  mean  speed  of  the 
ejected  particles  appears  to  saturate  at  roughly  10%  of  the  speed  of  the 
impacting  grain,  and  the  mean  ejection  angle  tends  to  be  oriented 
downwind  at  60-70°  from  the  horizontal.  All  results  for  impacts  between 
8°  and  15°,  which  cover  the  range  of  expected  impact  angles  in  eolian 
saltation,  display  only  slight  dependence  on  the  impact  angle.  These  results 
are  expressed  in  terms  of  the  splash  function  (Ungar  and  Haff,  1987),  an 
analytic  expression  for  the  number  density  of  grains  ejected  from  the  bed 
and  their  distribution  of  initial  velocities,  for  a  given  impact  velocity. 

Rebounds 

The  incident  grain  rebounds  with  probability  close  to  unity  for  the  impact 
speeds  simulated.  In  the  rare  cases  where  the  incident  grain  fails  to 
rebound,  the  micro-geometry  of  the  collision  pocket  is  such  that  the 
momentum  of  the  impact  is  efficiently  delivered  to  one  or  more  nearby 
grains,  while  the  colliding  particle  is  geometrically  trapped.  The  mean 
speed  of  the  rebound  is  60%  of  that  prior  to  collision  and  is  distributed  with 
a  standard  deviation  of  roughly  0.1  of  the  impact  speed.  The  ejection  angle 
is  significantly  steeper  than  the  impact  angle,  as  suggested  by  earlier  work 
of  Rumpel  (1985),  Werner  and  Haff  (1988),  and  others. 
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The  splash 


The  number  of  splashed  grains  resulting  from  an  impact  increases 
monotonically  with  the  impact  speed,  as  do  the  speed  and  angle  of  these 
ejections.  The  mean  number  of  ejections  increases  slightly  more  than 
linearly  with  impact  speed;  for  the  chosen  grain  characteristics,  impacts  of 
6m/s  result  in  10  ejecta.  The  mean  ejection  speed  for  the  splashed  grains 
increases  less  than  linearly  with  impact  speed,  becoming  approximately 
constant  (at  ~0.7m/s)  for  impacts  greater  than  5m/s.  The  mean  ejection 
angle  varies  little,  for  all  but  the  lowest  speed  impacts,  being  canted 
downwind  at  roughly  60-65°  form  the  horizontal.  Low  velocity  impacts 
result  in  either  no  ejecta,  or  ejecta  that  appear  very  similar  in  speed  and 
angle  with  that  of  a  rebounding  particle. 

Discussion  of  the  single  impact  simulations 

The  picture  of  the  grain  splash  emerging  from  the  numerical  experiments  is 
both  qualitatively  and  quantitatively  similar  to  that  arising  from  physical 
experiments  (e.g.,  Mitha  et  al.,  1986;  Willetts  and  Rice,  1986a,b,  1988),  and 
from  earlier  numerical  simulations  (e.g.,  Werner  and  Haff,  1986;  1988a).  A 
most  impressive  correspondence  comes  with  the  data  on  impacts  collected 
with  high  speed  motion  photography  taken  during  saltation  in  a  wind 
tunnel  (Willetts  and  Rice,  1986a,  1988). 

Although  considerable  caution  is  warranted  in  relying  upon  these 
numerical  experiments,  which  involve  two-dimensional,  perfect  circles,  we 
are  encouraged  by  the  correspondence  with  physical  experiments  using 
real  sand  grains  in  three-dimensional  pockets.  We  are  most  skeptical  of  the 
low  impact  velocity  results,  where  experimental  confirmation  is  difficult 
due  to  problems  associated  with  photography  very  near  the  bed.  At  these 
low  impact  speed  conditions,  the  three-dimensional  nature  of  the  real  beds 
should  lead  to  more  efficient  capture  of  low  velocity  grains,  leading  to  a 
reduction  in  the  probability  of  rebound. 

For  reasons  of  computational  efficiency,  the  simulated  beds  used  in  the 
impact  experiments  were  only  10  grain  diameters  wide;  they  tended 
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therefore  to  possess  microtopography  of  the  order  of  one  grain  diameter. 
This  contrasts  with  Bagnold’s  observation  that  saltation  impacts  create 
craters  with  rim  heights  on  the  order  of  two  grain  diameters.  In  addition, 
preliminary  measurements  of  microtopographic  profiles  with  lengths  of 
many  thousand  grain  diameters  on  actual  sand  beds  confirm  that,  even 
when  artificially  smoothed,  such  surfaces  have  micro-roughness  at 
wavelengths  much  greater  than  10D  with  amplitudes  of  up  to  4D. 

Impacts  with  the  upwind  sides  of  these  bumps  should  reduce  the  number 
of  rebounds  to  be  expected,  and  may  even  result  in  occasional  backwards- 
directed  ejecta  (or  rebounds)  that  we  will  have  missed  in  our  simulations, 
but  that  have  been  observed  in  physical  experiments  (B.B.  Willetts, 
personal  communication,  1988).  In  many  cases,  impacts  of  grains  with  the 
granular  surface  result  in  rebounds  that  appear  much  as  if  the  surface  was 
rigid  (Werner  and  Haff,  1988a),  the  angle  of  rebound  being  determined 
largely  by  the  surface  geometry,  specifically  the  orientation  of  the  surface 
normal.  A  knowledge  of  the  microtopographic  profiles  of  the  real  surfaces 
was  used  to  calculate  such  surface  normals  for  a  range  of  possible  impacts 
with  known  grain  size  and  impact  angle.  The  resulting  reflections  from  the 
surface  show  distributions  similar  to  those  resulting  from  our  calculated 
rebounds,  and  show  the  expected  dependence  of  the  distributions  on 
impact  angle  and  impacting  particle  diameter.  Although  our  calculations 
are  restricted  to  single  grain  size  beds  and  to  similarly  sized  impacting 
grains,  these  profiles  allow  estimation  of  the  effects  of  variable  impacting 
particle  diameters.  Larger  grains  impacting  a  finer  surface  will  perceive  a 
smoother  surface,  and  should  result  in  a  lower  angle  rebound,  while  finer 
grains  encountering  a  coarser  surface  ought  to  rebound  at  greater  angles. 
The  geometrical  argument  presented  here  should  be  altered  by  the  fact 
that  larger  particles  will  create  larger  disturbances  of  the  surface; 
simulations  of  these  events  are  needed  to  clarify  the  relative  roles  of  the 
geometrical  and  dynamical  effects. 

The  statistics  of  the  rebounds  and  splashed  grains  are  effectively  mingled 
for  low  impact  velocity  events,  and  become  more  distinguishable  at  higher 
velocities.  In  all  cases,  however,  an  event  may  be  most  easily  characterized 
by  the  sum  of  these  two  separate  populations,  one  representing  the 


rebound,  the  other  the  splashed  grains.  We  follow  this  approach  in 
condensing  the  results  of  these  simulations  into  a  useful  analytical 
representation  of  the  grain-impact  process  to  be  used  in  the  saltation 
simulations  to  follow. 

Splash  function 

The  results  of  the  single  impact  simulations  are  condensed  into  analytical 
expressions  for  the  number  of  ejecta  and  the  probability  distribution  of 
their  ejection  velocities,  as  functions  of  impact  velocity  (Anderson  and  Haff, 
1988).  As  we  found  little  dependence  of  the  properties  of  the  grain  splash 
over  the  relevant  range  of  expected  angles  of  impact  in  eolian  saltation,  we 
have  ignored  any  dependence  on  impact  angle.  A  reasonable  expression  for 
the  number  of  rebounding  particles  to  be  expected  from  a  single  impact  of 
velocity  Vjmi  in  each  of  j  ejection  velocity  "bins",  labelled  VQj,  is  a  gaussian: 

Nr(Voj)  =  a  exp  (-(Voi-bVimi)/cVimi)2)  dVo 

where  a,  b  and  c  are  constants. 

The  distribution  for  the  expected  splashed  grains  is  an  exponential,  this 
time  scaled  by  the  expected  number  of  grains  ejected  per  impact,  fVimi 

Ne(Voj)  =  [fVimi]  exp  (-Voi/hVim)  dVo 

In  each  case  dVo  is  the  width  of  the  ejection  velocity  bin.  The  total  number 
of  grains  leaving  the  bed  in  each  velocity  bin  subsequent  to  a  single  impact 
of  velocity  Vimi,  then,  is  the  sum  of  these  two  expressions:  Ns  =  Nr+Ne. 

Modification  of  the  wind  velocity  profile 

As  grains  are  accelerated  by  the  force  of  the  wind,  they  impose  an  equal 
and  opposite  force  on  the  wind.  The  profile  of  this  force  per  unit  volume  on 
the  wind  was  calculated  and  used  to  alter  the  effective  stress  available  to 
shear  the  air  at  all  levels  in  the  flow.  Prediction  of  particle  concentration, 
mass  flux,  and  kinetic  energy  flux  caused  by  wind  requires  accurate 
calculation  of  particle  trajectories,  which  in  turn  requires  a  detailed 
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knowledge  of  the  wind  velocity  profile.  A  theoretical  model  for  the 
magnitude  and  distribution  of  the  drag  imposed  upon  the  wind  by  saltating 
grains  is  necessary  to  establish  the  relative  magnitude  of  the  three 
contributions  to  the  overall  drag  over  a  mobile  sand  bed:  saltating  grains, 
stationary  grains  on  the  bed,  and  form  drag  due  to  ripples.  In  the  eolian 
sediment  transport  system,  the  vertical  region  within  which  the  velocity 
profile  is  modified  by  sediment  transport  is  of  the  order  of  centimeters, 
making  detailed  measurement  of  the  velocity  structure  possible. 

In  the  lowest  10  m  of  the  atmospheric  boundary  layer  the  shear  stress,  Xb, 
is  approximately  constant.  Within  this  region,  but  well  above  the  heights  of 

the  roughness  elements,  a  logarithmic  velocity  profile  is  expected  to  exist, 
dependent  upon  a  single  velocity  scale,  u*,  and  a  single  length  scale,  z0.  The 

velocity  scale  is  the  shear  velocity,  u*,  defined  as  Vtb/pb,  where  pa  is  the  air 

density.  The  shear  stress,  in  turn,  is  governed  by  larger  scale  atmospheric 
circulation,  driven  by  pressure  gradients  imposed  largely  by  differential 
solar  heating  of  the  atmosphere.  In  the  absence  of  large  bedforms,  the 
length  scale  setting  the  logarithmic  velocity  profile  is  proportional  to  the 
diameter  of  the  grains  in  the  bed:  z0  =  D/30  for  uniform  grains. 

Sediment  transport  alters  the  effective  roughness  of  the  bed  in  two  ways: 
first,  the  horizontal  acceleration  of  transported  grains  extracts  momentum 
from  the  wind,  and  second,  the  formation  of  ripples  in  the  bed  imposes  a 
form  drag  on  the  wind.  The  resulting  sediment  transport  roughness,  zost, 

measured  by  extrapolating  the  velocity  profile  outside  the  sediment 
transport  region  to  the  U=0  axis,  has  been  shown  to  depend  on  the  shear 
velocity  [after  Owen  (1964)] 


u*2 

Z°st~a^i_ 


with  a~  .02. 

Because  we  expect  the  mean  vertical  velocity  of  particles  upon  ejection 
from  the  bed  to  increase  with  the  shear  velocity,  and,  as  argued  above,  in 
the  absence  of  non-gravitational  forces,  a  grain  leaving  the  bed  with 
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velocity  u*  reaches  a  height  of  u*2/2g,  the  success  of  Owen’s  formulation 

implies  that  the  roughness  is  no  longer  proportional  to  the  size  of  the 
stationary  roughness  elements  in  the  bed,  but  rather  to  some  typical  hop 
height  of  a  saltating  grain. 

Wind  velocity  measurements  within  the  saltation  region  show  significant 
deviations  from  the  logarithmic  profile.  Velocity  gradients  nearest  the  bed 
are  reduced  from  those  expected  by  extrapolating  the  outer  flow  toward 
the  bed,  and  result  in  ln(z)  -  U  plots  of  wind  profiles  that  are  convex 
upward.  This  curvature  corresponds  to  Bagnold's  (1941)  "kink"  in  the 
profile,  which  occurs  typically  on  the  order  of  one  to  a  few  centimeters 
above  the  bed.  Such  near-bed  deviations  from  the  logarithmic  velocity 
profile  during  sediment  transport  have  led  to  incorrect  assessment  of  the 
shear  velocity  during  experimental  work  in  wind  tunnels.  All  too  often,  the 
slope  of  the  least-squares  fit  to  the  ln(z)-U  plot,  from  which  the  shear 
velocity  is  calculated,  incorporates  at  least  in  part  the  region  in  which 
saltating  grains  are  expected  to  impose  a  systematic  departure  from  the 
simple  logarithmic  profile.  Correct  assessment  of  the  shear  velocity  is 
essential,  for  instance,  in  the  development  of  the  correct  functional 
dependence  of  the  total  mass  flux  on  the  shear  velocity  (e.g.,  Bagnold,  1941 
and  compiled  in  Greeley  and  Iversen,  1985). 

A  full  model  of  the  modification  of  the  wind  profile  by  saltating  grains 
must  therefore  be  able  to  predict  both  this  curvature  of  the  profile,  and 
the  altered  effective  roughness  of  the  bed.  As  emphasized  by  Ungar  and 
Haff  (1987),  a  complete  steady  state  saltation  model  must  include  iteration 
through  a  wind-velocity  profile  feedback  loop,  as  the  altered  wind  profile 
will  change  the  suite  of  particle  trajectories,  which  in  turn  result  in  a  new 
wind  profile,  etc. 

Although  it  has  long  been  recognized  that  the  modification  of  the  wind 
profile  results  from  the  extraction  of  momentum  from  the  wind  by  the 
saltating  grains,  previous  attempts  to  incorporate  this  effect  have  suffered 
from  oversimplification  of  the  saltation  process.  Most  notably,  Owen  (1964) 
assumed  that  all  particles  trace  identical  trajectories.  It  is  now  clear. 
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however,  that  the  probabilistic  nature  of  the  grain-bed  interaction  leads  to 
a  broad  distribution  of  initial  conditions  for  particle  trajectories,  and  to 
large  gradients  in  particle  concentration,  mass  flux,  and  kinetic  energy  flux 
with  height  above  the  bed  (e.g.,  Anderson  and  Hallet,  1986;  Anderson, 
1986a, b;  Jensen  and  Sorenson,  1985;  Sorenson,  1986;  Willetts  and  Rice, 
1986a,b).  This  structure  is  expected  to  be  reflected  in  the  profile  of  the 
force  imposed  on  the  wind  by  the  horizontal  acceleration  of  transported 
grains. 

Two  recent  attempts  have  been  made  to  calculate  wind  profiles  during 
sediment  transport  (Ungar  and  Haff,  1987;  Sorensen,  1986).  Ungar  and  Haff, 
in  their  analysis  of  the  saltation  problem,  confine  their  calculations  to  the 
simplest  possible  case  that  retains  all  the  important  elements  of  the 
problem.  At  any  shear  velocity,  they  force  their  solution  to  yield  only  one 
particle  trajectory.  As  the  grain-bed  interaction,  characterized  by  their 
"splash  function",  is  independent  of  wind  speed,  i.e.  the  same  liftoff  velocity 
is  retained  for  a  particular  grain  impact  velocity  no  matter  what  the  wind 
structure  is,  the  single  trajectory  allowed  at  each  shear  velocity  must  have 
the  same  impact  velocity.  This  requires  that  each  such  trajectory 
experience  the  same  net  acceleration  by  the  wind.  They  argue,  therefore, 
that  in  accord  with  Bagnold  (1941)  "no  matter  how  hard  the  wind  is  made 
to  blow  ...  the  wind  velocity  at  a  height  of  about  3  cm  remains  almost  the 
same.  Moreover,  at  levels  still  closer  to  the  ground  the  wind  velocity 
actually  falls  as  the  wind  above  is  made  stronger."  Ungar  and  Haff’s 
computed  profiles  show  this  behavior.  It  remains  to  be  seen  whether 
similar  results  can  be  obtained  for  a  more  realistic  range  of  particle 
trajectories,  arising  from  a  more  realistic  grain-bed  interaction,  or  "splash 
function." 

Sorenson's  treatment  of  wind  velocities  during  sediment  transport  allows  a 
realistic  range  of  particle  trajectories.  The  principal  contrast  with  the 
formalism  presented  below  lies  in  the  nature  of  the  postulated  "closure" 
relation  between  the  fluid  stress  and  the  shear  rate,  discussed  further 
below. 
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Momentum  equation  for  the  air 


We  seek  an  expression  for  the  momentum  of  the  air  during  saltation  that 
both  takes  account  of  the  momentum  extracted  by  the  grains  during 
saltation,  and  reduces  properly  to  the  "law  of  the  wall"  appropriate  for  the 
planetary  boundary  layer  in  the  absence  of  saltation.  Following  the  work 
of  Ungar  and  Haff,  we  first  formalize  the  effect  of  saltating  grains  on  the  air, 
which  results  in  the  establishment  of  an  extra  body  force  term  in  the 
Navier-Stokes  equation  for  the  force  balance  on  a  parcel  of  air.  We  then 
seek  an  appropriate  form  for  the  distribution  of  this  body  force  with  height 
above  the  bed.  Finally,  we  suggest  a  possible  closure  relating  the  turbulent 
stresses  to  the  local  mean  wind  shear  that  remains  consistent  with  the  law 
of  the  wall  in  the  absence  of  saltation. 

Within  the  saltating  curtain,  the  assumption  of  constant  shear  stress  used  to 
construct  the  law  of  the  wall  breaks  down,  as  the  acceleration  of  massive 
grains  imposes  an  additional  force  on  the  wind.  Following  the  formulation 
of  Ungar  and  Haff,  we  identify  a  horizontal  body  force  on  the  wind,  Fx(z), 

acting  in  the  upwind  (negative-x  )  direction  due  to  the  acceleration  of  the 
grains  by  the  air.  This  appears  explicitly  in  the  turbulent  Navier-Stokes 
equations  as  an  additional  force  term  operating  in  the  negative-x  direction: 

3u  ~ 

Pa-^r  +  PaU*  Vu  =  -  Vp  +  V-xt  -  pa  g  -  Fx 

where  pa  is  the  air  density,  u  is  the  mean  horizontal  wind  velocity,  g  is  the 
acceleration  due  to  gravity,  and  xt  is  the  turbulent  (Reynolds)  shear  stress, 
representing  the  flux  of  momentum  via  correlations  in  the  velocity 
fluctuations.  Given  steady  (9/9t~0),  horizontally  uniform  flow  (u-Vu=0),  and 
making  boundary  layer  approximations  (9/9z  »  9/9x,  9/9y),  the  equation  for 
momentum  in  the  downwind  (positive-x  )  direction  collapses  to: 

=  Fx(z) 
dz 
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In  the  absence  of  saltating  grains  the  right  hand  side  vanishes,  and  the  first 
integration  yields  xt  =  constant.  After  identifying  this  constant  as  the  shear 

stress  imposed  by  the  exterior  flow,  x\y  =  pa  11*2,  making  the  common 

"closure  hypothesis"  that  the  turbulent  stresses  may  be  identified  as  the 
product  of  an  eddy  diffusivity,  K,  with  the  strain  rate,  9u/8z,  and  making 

the  further  assumption  that  the  eddy  diffusivity  varies  linearly  with 
height,  K=ku*z,  where  k  is  von  Karman’s  constant  (=0.40),  a  second 

integration  yields  the  well  known  law  of  the  wall,  or  logarithmic  profile: 
u  =  (u*/k)ln(z/z0),  where  z0  is  the  effective  roughness  of  the  bed.  Such 

conditions  should  apply  throughout  the  profile  in  the  absence  of  sediment 
transport,  and  above  the  region  within  which  grains  are  being  appreciably 
accelerated  by  the  wind  during  sediment  transport,  the  difference  between 
them  being  in  the  value  of  the  effective  roughness,  z0. 

Within  the  saltation  region,  however,  the  stress  on  the  wind  must  vary  in 
the  vertical  direction  as  the  force  on  the  wind  due  to  the  extraction  of 
momentum  by  saltating  grains  varies.  Assuming  a  constant  total  stress 
available  for  transporting  momentum  of  either  grains  or  fluid  across  any 
level  z,  the  stresses  may  be  partitioned  according  to 

IZmax 

Fx(z)  dz 

where  zmax  is  the  maximum  height  to  which  a  saltating  grain  travels  in  the 
given  transport  conditions. 

The  first  term  on  the  right  hand  side  represents  the  stress  available  to 
shear  the  air  at  the  level  z  ,  or  the  flux  of  fluid  momentum  across  that  level; 
the  second  ,  also  called  the  "grain  stress",  xg  (Sorenson,  1986),  represents 
the  change  in  horizontal  momentum  of  grains  between  their  upward  and 
downward  crossings  of  the  level  z  ,  or  the  net  flux  of  horizontal  grain 
momentum  across  that  level.  We  see  that  as  the  bed  is  approached  from 
above,  the  grain  stress  increases,  leaving  less  shear  stress  available  to  shear 
the  fluid. 
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This  equation  may  be  rewritten  as 


tt(z) 


|  Fx(z)dz  +Tsf 

0 


where  xsf  is  the  fluid  stress  at  the  bed,  or  the  skin  friction,  part  of  which 
may  be  taken  up  by  form  drag  due  to  ripples. 


We  see  then  that  the  skin  friction,  or  the  shear  stress  exerted  on  the  bed  by 
the  wind,  is  simply  the  far-field  shear  stress,  Xb,  minus  the  total  change  in 
horizontal  momentum  of  all  grains  ejected  from  a  unit  area  of  bed  in  a  unit 
time 


*sf  — 


Fx(z)dz 


For  a  given  exterior  wind  condition,  characterized  by  xb,  as  the  body  force 
diminishes  in  magnitude,  more  fluid  stress  is  made  available  at  the  bed  for 
aerodynamic  initiation  of  saltation;  and  conversely,  as  body  force  increases, 
the  shear  at  the  bed  decreases.  Owen  (1964)  hypothesized  that  the  self- 
regulatory  nature  of  eolian  saltation  arose  from  such  a  feedback,  with  skin 
friction  held  near  the  threshold  shear  stress  necessary  to  entrain  particles 
aerodynamically.  It  remains  to  test  this  hypothesis  quantitatively. 


We  now  seek  an  equation  for  the  velocity  gradient  as  a  function  of  height, 
which  when  integrated  will  yield  a  velocity  profile  in  equilibrium  with 
sediment  tranport.  Given  the  above  relations  between  far  field  shear 
stress,  skin  friction,  and  grain  stress,  we  require  both  a  constitutive  relation 
between  the  turbulent  stresses  and  the  mean  velocity  gradient,  and  the 
body  force  profile,  Fx(z).  For  simplicity,  we  again  postulate  an  eddy 
diffusivity  closure:  xt=paK9u/3z,  and  retain  the  linear  dependence  of  the 
eddy  viscosity  with  height,  K=ku*z  ,  that  gives  rise  to  the  logarithmic 

velocity  profile  in  the  absence  of  sediment  transport.  However,  referencing 
u*  to  the  total  stress,  x^,  is  no  longer  appropriate.  We  define  a  local ,  or 
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effective  shear  velocity,  referenced  to  the  local  turbulent  stress,  xt,  as 
U*eff(z)=(xt/Pa Then,  from  xt  =  pakzu*eff9u/9z,  we  have 


%  -  Fx(z)dz 

du  _  i 

L 

dz  kz 

L  Pa  J 

Note  that  above  the  region  within  which  grains  accelerate,  i.e  when  z  >  zmax, 
or  in  the  absence  of  sediment  transport  altogether,  i.e.  when  Fx=0  for  all  z, 
the  numerator  becomes  (xb)1/2,  and  the  rate  of  shear  becomes  du/dz  = 
u*/kz,  which  again  yields  the  logarithmic  profile,  as  required.  Given  the 

form  of  the  force  profile  at  any  time  in  the  evolution  of  the  saltation 
population,  derived  in  the  next  section,  this  equation  may  be  numerically 
integrated  to  yield  the  wind  velocity  profile  throughout  the  saltation  region. 

Note  that  irrespective  of  the  actual  form  of  the  force  profile,  the  fluid 
stress,  xt,  and  hence  the  effective  shear  velocity,  will  increase  monotonically 
with  height,  giving  rise  to  convex  upward  u-ln(z)  plots.  Those  profiles 
showing  inflections  in  semi-log  space,  as  some  of  Bagnold's  (1941)  early 
profiles  do,  are  suspect;  the  stress  available  to  shear  the  fluid,  and  hence 
the  effective  shear  velocity,  should  increase  monotonically  away  from  the 
bed,  yielding  convex-upward  wind  velocity  profiles  in  semi-log  space. 
Inflections  resulting  from  an  initial  decrease  in  effective  shear  velocity, 
followed  by  an  increase  in  shear  velocity,  should  not  exist. 

Force  on  the  wind  due  to  saltating  particles 

Following  the  formalism  established  in  previous  work  (Anderson  and 
Hallet,  1986;  Anderson,  1986),  the  force  due  to  identical  trajectories  with 
unit  ejection  rate  is  first  calculated,  i.e.,  we  compute  numerically  the 
ejection  Greens  function;  the  distribution  of  initial  conditions  and  the  actual 
ejection  rate  at  any  time  in  the  evolution  of  the  saltation  population  are 
then  incorporated  to  yield  a  total  instantaneous  force  profile. 
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The  highest  instantaneous  force  is  attained  early  in  the  ascending  limb  of 
the  trajectory,  before  the  particle  has  been  appreciably  accelerated  by  the 
wind,  and  where  the  relative  velocity  between  the  particle  and  the  air  is 
therefore  the  greatest.  The  force  on  the  particle  becomes  negative  shortly 
before  impacting  the  bed,  indicating  that  the  particle  there  is  travelling 
faster  than  the  wind. 

Summing  over  the  ascending  and  descending  limbs  of  the  trajectory,  and 
assuming  a  single  particle  is  ejected  per  unit  area  of  bed  per  unit  time,  i.e. 
Nj=l,  with  initial  vertical  velocity  w0,  results  in  the  "identical  trajectory" 

force  on  the  wind: 


Fx(zlw0,Ni)  =  NiM 


[  t(w(w0))|+  |(w(w0))|.J 


where  <w>  is  the  mean  vertical  particle  velocity  in  crossing  the  height 
element  (z-dz,  z),  and  the  +  and  -  denote  upward  and  downward  limbs  of 
the  trajectory,  respectively  (Ungar  and  Haff,  1987;  Anderson,  1986b).  The 
upwind  direction  of  the  force  is  left  implicit.  The  "identical  trajectory"  force 
profile  displays  a  distinct  maximum  at  the  top  of  the  trajectory,  similar  to 
the  maxima  in  particle  concentration,  mass  flux,  and  kinetic  energy  flux 
profiles  reported  earlier  (Anderson  and  Hallet,  1986;  Anderson,  1986a). 
Although  the  instantaneous  horizontal  force  on  the  particle  peaks 
approximately  one  third  of  the  way  up  the  ascending  limb  of  the  trajectory, 
most  of  the  particle's  horizontal  acceleration  occurs  near  the  top  of  the  hop, 
where  its  vertical  velocity  is  very  low,  and  therefore  where  it  spends  the 
most  time.  A  similar  argument  explains  the  maximum  at  the  top  of  the 
particle  path  in  the  concentration,  mass  flux,  and  kinetic  energy  flux 
profiles  for  "identical  trajectory"  models  (Anderson  and  Hallet,  1986; 
Anderson,  1987). 

The  probability  density  of  the  vertical  liftoff  velocity,  p(w0)dw0,  and  the 
actual  number  of  particles  ejected  per  unit  area  of  bed  per  unit  time,  N,  are 
now  introduced  to  yield  an  integral  equation  for  the  total  horizontal  body 
force  per  unit  volume  on  the  wind  due  to  the  presence  of  saltating  particles: 
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f*(z)=rt1 


p(wo)  Fx(zlw0,Ni)  dw0 


In  earlier  work  (Anderson,  1986b),  force  profiles  and  the  resulting  wind 
velocity  profiles  were  calculated  using  reported  total  mass  fluxes  and 
assumed  forms  for  the  probability  distribution,  p(wQ),  to  constrain  the 

values  of  the  total  ejection  rate,  N,  in  the  above  equation.  These  calculations 
revealed  that  the  effect  of  the  form  of  the  probability  distribution  of  liftoff 
velocities  dominates  over  the  shape  of  the  "identical  trajectory"  force 
profile  in  determining  the  result  of  the  integration.  Force  profiles  decline 
sharply  above  the  bed,  with  a  scale  height  that  is  on  the  order  of  l-3cm. 

The  effect  of  transported  grains  on  the  wind  profile,  for  moderate  winds, 

should  essentially  vanish  by  5cm  above  the  bed,  as  is  observed.  As  scale 
height  is  a  reflection  largely  of  the  probability  distribution  of  initial 
velocities,  one  expects  this  will  be  strongly  influenced  by  grain  size,  and  by 
the  external  forcing  characterized  by  u*. 

The  choice  of  a  linearly  increasing  eddy  diffusivity  within  the  saltation 
region  follows  more  closely  the  treatments  of  Ungar  and  Haff  (1988)  and  of 
Sorenson  (1986),  (and  is  identical  to  the  mixing  length  algorithm  used  by 
Werner,  1990)  and  diverges  from  that  of  Owen  (1964),  who  argued  that  a 
constant  eddy  diffusivity  would  be  appropriate.  Owen  made  the  plausibility 
argument  that  the  intensity  of  turbulence  and  the  mixing  lengths  of  the 
turbulence  ought  both  to  be  dominated  by  wakes  cast  by  saltating 
particles,  and  should  therefore  be  roughly  constant  within  the  saltation 

curtain.  However,  the  concentration  of  particles  in  the  flow  is  on  the  order 
of  10-2  to  10*4  near  the  bed  (Anderson  and  Hallet,  1986;  Sorenson,  1986). 
It  may  therefore  be  expected  that  the  nearby  presence  of  a  continuous 
rough  bed  will  dominate  over  the  wakes  cast  by  these  sparsely  distributed 
particles  in  setting  the  length  scale  of  the  turbulence.  The  use  of  a  linearly 
varying  eddy  diffusivity  also  allows  a  much  simpler  feathering  of  the 

saltation  region  with  the  constant  stress  region  above,  rather  than  requiring 
the  abrupt  change  in  the  nature  of  the  turbulence  suggested  by  Owen.  This 
explicitly  recognizes  that  the  saltation  region  is  neither  physically  isolated. 
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nor  even  easily  identified  in  the  real  world.  The  character  of  the  saltation 
curtain  is  such  that  all  profile  quantities,  including  the  force  on  the  wind 
imposed  by  the  accelerating  particles,  fall  off  rapidly  away  from  the  bed,  as 
a  consequence  primarily  of  the  heavily  skewed  nature  of  the  probability 
density  of  liftoff  conditions  resulting  from  the  complex  grain-bed 
interaction.  It  is  a  sediment  transport  boundary  layer. 

That  the  shear  velocity  be  tied  to  the  local  fluid  stress,  x t,  rather  than  the 
total  or  far-field  stress,  X5,  makes  the  present  treatment  different  from 
those  of  both  Ungar  and  Haff  (1987),  and  Sorenson  (1986).  Such  an 
hypothesis  has  proven  fruitful  in  the  analysis  of  aqueous  systems,  where 
the  form  drag  due  to  multiple  sets  of  bedforms  is  modelled  with  a 
corresponding  number  of  matched  logarithmic  profiles,  each  characterized 
by  a  shear  velocity  referenced  to  the  spatially  averaged  form  drag 
associated  with  that  particular  bedform  scale.  The  present  expression  for 
the  wind  shear  identifies  both  the  magnitude  of  the  stress  due  to  saltating 
grains,  and  its  profile,  allowing  a  smoothly  varying  wind  profile  rather  than 
a  matched  logarithmic  profile. 

A  further  difficulty  with  Owen's  formulation  of  the  wind  profile  within  the 

saltation  region  is  that  the  no-slip  condition  is  left  unsatisfied.  Using  his 
formula,  wind  velocities  at  the  bed  remain  of  the  order  7-8  u*.  Although 

this  has  prompted  earlier  workers  to  justify  use  of  a  constant  wind  velocity 
(=8.5  u*)  in  the  saltation  layer,  thereby  simplifying  trajectory  calculations, 

the  present  formulation  is  considerably  more  realistic,  is  not 
computationally  taxing,  and  provides  a  profile  valid  through  the  entire 
region  of  interest  in  sediment  transport  mechanics. 

Saltation  model 

The  actual  particle  motion  through  the  air  is  handled  as  follows.  Because  of 
computational  limitation,  particles  are  allowed  to  leave  the  bed  with  one  of 
only  ten  distinct  velocities;  i.e.,  the  distribution  of  initial  conditions  is 
discretized.  These  are  chosen  such  that  the  lowest  liftoff  velocity  trajectory 
reaches  a  height  8  equal  to  approximately  one  grain  diameter,  while  that 
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with  the  highest  liftoff  velocity  reaches  a  height  above  which  little  saltation 
flux  is  expected.  This  latter  limit  was  chosen  to  be  wo=5u*,  meaning  that 

typical  trajectory  heights  for  these  grains  is  8=25u*2/2g,  or  roughly  1.2m 
for  u*  =  1.0m/s.  The  intervening  liftoff  velocity  "bins"  are  distributed 

logarithmically  in  order  to  maximize  the  profile  information  in  the  lower 
portion  of  the  saltation  curtain,  where  most  of  the  grains  travel.  At  any 
instant  in  time,  each  liftoff  velocity  bin  (V0)j,  contains  Nj  particles. 


Collision  list 

Participation  in  the  splash  process  requires  impact  with  the  surface,  which 
occurs  with  decreasing  frequency  as  the  initial  velocity  of  the  trajectory 
increases.  (As  noted  previously,  trajectory  times  are  well  estimated  by 
thop=b(2w0/g),  where  wo  is  the  vertical  component  of  the  liftoff  velocity, 

and  b  is  a  coefficient  of  value  close  to  1,  only  weakly  dependent  upon  the 
initial  velocity  of  the  grain.)  To  account  for  the  variable  trajectory 
durations,  a  "collision  list"  is  established.  For  each  trajectory  class  a  list  is 
maintained  of  the  number  of  grains  expected  to  be  impacting  the  surface  at 
each  time  increment  in  the  future.  At  each  time  step,  then,  the  collision  list 
for  each  trajectory  class  is  interrogated  to  reveal  the  number  of  impacts  in 
each  trajectory  bin.  Their  impact  velocities  are  then  assessed  by  calculating 
the  trajectories  using  the  (changing)  wind  velocity  profile  averaged  over 
the  hop  time  of  the  trajectory;  the  number  of  splashed  particles  in  each  of 
the  liftoff  velocity  bins  is  then  calculated  from  the  splash  function.  Finally, 
the  collision  list  is  updated  in  accord  with  this  new  set  of  ejecta,  and  the 
entire  collision  list  is  moved  forward  in  time  one  step  to  reveal  the  new 
numbers  of  impacts  in  each  bin  to  be  expected  upon  the  next  time  step. 

In  the  calculation,  time  is  discretized  at  two  levels.  The  time  step 
associated  with  the  reassessments  of  the  wind  profile  and  the  grain  splash 
is  determined  by  the  shortest  trajectory  time,  or  roughly  2(2D/g)1/2.  A 
second,  much  shorter  time  step  is  needed  to  assure  that  each  trajectory  is 
calculated  with  sufficient  precision  to  account  correctly  for  the  impact 
velocity  and  the  various  profile  quantities  of  concern.  Were  we  to  have  to 
calculate  the  actual  grain  dynamics  of  the  bed  upon  each  and  every  impact 
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with  the  bed,  a  third,  much  shorter  time  step  would  be  involved 
corresponding  to  the  elastic  material  properties  of  the  bed  grains;  this 
would  be  of  the  order  of  10'^  seconds.  The  use  of  the  splash  function  to 
represent  the  response  of  the  bed  to  impacts  removes  this  necessity. 

Full  simulations 

Inputs  required  for  the  full  simulation  of  the  eolian  saltation  system  are  (i) 
a  grain  diameter  and  density,  (ii)  an  initial  wind  velocity  profile,  the 
constants  determining  the  splash  function,  and  the  coefficient  of 
proportionality  between  the  excess  shear  stress  and  the  aerodynamic 
entrainment  rate,  a.  The  initial  wind  velocity  profile  is  taken  to  be 
logarithmic.  For  a  flat  bed,  the  roughness  height  is  related  to  the  size  and 
packing  of  the  stationary  sand  grains  in  the  bed  (=D/30).  In  his 
experiments  on  the  effect  of  sand  movement  on  the  surface  wind,  Bagnold 
(1941)  produced  wind  profiles  over  a  wetted  sand  bed  that  was  previously 
"not  only  pitted  with  tiny  bombardment  craters  a  few  grain  diameters  in 

size,  but  was  made  to  undulate  in  the  usual  flat  transverse  ripples".  The 
resulting  u-  ln(z)  profiles  for  u*=.20-.62m/s  show  little  or  no  curvature 

between  2mm  and  10cm  above  the  bed,  and  all  yield  roughness  heights 
closely  approximated  by  zo=2D/30,  where  2D  is  roughly  the  mean  height  of 

the  impact  craters  (Bagnold,  1941).  Accordingly,  the  initial  roughness 
height  was  chosen  to  be  2D/30.  The  shear  velocity,  u*,  in  turn  sets  the 

initial  shear  stress  at  the  bed,  (xa)0=Tb=pau*2,  to  be  used  in  subsequent 
recalculations  of  the  stress  and  wind  profiles.  Simulations  are  run  until  a 
steady  state  is  achieved,  a  state  characterized  by  little  or  no  change  in  the 
mass  flux,  the  wind  velocity  profile,  or  in  the  collision  list  from  one  time 
step  to  the  next. 

In  a  typical  simulated  evolution  of  the  saltation  population  and  of 
associated  quantities,  initial  entrainment  is  entirely  aerodynamic.  These 
grains  gain  horizontal  momentum  from  the  wind,  and  impact  with  velocities 
such  that  a  small  proportion  rebound  with  greater  initial  velocities  than  the 
initial  trajectories;  few  subsequent  grains  are  splashed  at  this  stage.  An 
initial  delta  function  probability  distribution  of  liftoff  velocities  therefore 
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evolves  to  a  broader  probability  distribution,  filling  out  into  the  higher 
velocities.  The  grains  with  higher  liftoff  velocities  are  airborne  for  longer 

periods  of  time  before  impacting  and  contribute  more  strongly  to  the  splash 

ejection  of  other  grains.  The  full  range  of  possible  ejection  velocities  is 
populated  only  after  many  tens  of  short-trajectory  times.  At  this  point  the 
total  number  of  grains  in  transport  begins  to  grow  rapidly,  the  highest 
impact  velocity  grains  being  most  efficient  in  splashing  grains  into  the 
airstream.  The  resulting  roughly  exponential  growth  is  curtailed  only  when 

the  extraction  of  momentum  from  the  wind  is  sufficient  to  alter 

significantly  the  wind  velocity  profile,  which  in  turn  alters  the  impact 
velocities  of  the  longer  grain  trajectories. 

The  system  eventually  reaches  a  steady  state  characterized  by  a  specific 
total  mass  flux,  an  equal  number  of  impacting  and  ejected  grains,  and  a 
stationary  wind  velocity  profile.  The  overshoot  of  the  steady  state  appears 
to  be  due  to  the  time  lag  associated  with  the  0.2  to  0.3  second  hop  times  of 
the  most  energetic  trajectories  responsible  for  the  majority  of  the  ejections 
from  the  bed.  The  steady  state  mass  flux  is  well  within  the  range  of  fluxes 
measured  in  wind  tunnel  studies  for  the  same  combination  of  shear 
velocity  and  grain  size.  Whereas  the  initial  saltation  population  is  entirely 
aerodynamic,  the  steady  state  saltation  population  for  most  imposed  shear 
velocities  contains  no  aerodynamically  entrained  grains.  As  the  population 
of  splashed  and  rebounding  grains  increases,  the  shear  stress  at  the  bed  is 
reduced,  which  produces  a  corresponding  decrease  in  the  rate  of 
aerodynamic  entrainment.  The  steady  state  shear  stress  at  the  bed  is 
reduced  to  slightly  below  the  critical  shear  stress  for  aerodynamic 
entrainment.  The  effective  roughness  of  the  bed  is  simultaneously  greatly 
increased,  as  reflected  in  the  rise  of  the  U=0  intercept  (from  z0  to  z0'  ).  This 

change  in  bed  roughness  is  in  accord  with  numerous  measured  wind 
velocity  profiles  during  saltation  experiments,  as  summarized  by  Owen.  The 
resulting  steady  state  profile  of  mass  flux  displays  the  characteristic  rapid 
decrease  above  the  bed,  implying  that  the  system  has  evolved  to  produce  a 
realistic  probability  distribution  for  the  initial  trajectory  velocities. 

By  varying  the  initial  shear  velocity,  several  such  simulations  were 
performed  (Anderson  and  Haff,  1988)  to  produce  a  mass  flux  relation  that 
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is  broadly  similar  to  those  derived  empirically:  above  a  threshold  shear 
velocity  (u*c=  (tc/pa)1/2),  the  flux  increases  as  a  power  of  (u*-u*c). 

Further  calculations  demonstrated  yet  another  long-recognized  feature  of 
eolian  saltation:  hysteresis.  Once  steady  state  had  been  established  using  a 
particular  shear  velocity,  the  externally  imposed  shear  stress  (and  hence 
u*)  was  diminished  (in  stepped  changes)  to  below  that  necessary  to  entrain 

grains  aerodynamically.  For  choices  of  this  shear  velocity  down  to 
approximately  0.75  of  the  aerodynamic  threshold  shear  velocity,  u*c,  a  low 

but  steady  and  finite  flux  was  achieved.  This  reflects  the  hysteresis  in  the 
system  that  forced  Bagnold  to  define  two  threshold  velocities:  a  fluid 
threshold,  and  an  impact  threshold,  the  latter  being  on  the  order  of  0.8  of 
the  former. 

The  response  time  of  the  saltation  system  appears  to  be  on  the  order  of  1  to 
2  seconds,  or  several  long-trajectory  times  (Anderson  and  Haff,  1988; 
Anderson,  1988).  For  the  cases  run  to  date,  it  appears  that  the  response 
time  is  a  weak  function  of  the  shear  velocity,  being  longer  for  lower  shear 
velocities.  A  knowledge  of  the  response  time,  which  is  difficult  to  obtain 
from  wind  tunnel  experiments,  will  allow  us  to  treat  the  problem  of 
predicting  total  mass  fluxes  in  variable  winds. 
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